# -*- coding: utf-8 -*-

# %%
import xlrd
import numpy as np
from matplotlib import pyplot as plt
import os
import time
import itertools
import collections
import Quench_Initial_Settings_zdir
import sys
from scipy import ndimage as nd

sys.path.append(r'X:\Work\Rb87Python')
sys.path.append(r'D:\vc\Rb87Python')

import basic_function.Rb87BEC_Basic_Functions as bf
import basic_function.Rb87BEC_Fit as bffit

import socket

os.system('chcp 65001')
Host_Name = socket.gethostname()
if Host_Name == 'VictorChengPC':
    Computer_Name = 'X:\\Work\\mail.ustc.edu.cn\\Rb87BEC - Files\\'
elif Host_Name == 'VictorChengXPS':
    Computer_Name = 'X:\\Work\\mail.ustc.edu.cn\\Rb87BEC - Files\\'
elif Host_Name == 'calcenter2':
    Computer_Name = 'D:\\vc\\mail.ustc.edu.cn\\Rb87BEC - Files\\'

# Index = 0.05

# %%
time_start = time.time()
time_start_struct = time.localtime(time_start)

# %% Switches
FitBEC_Switch = 0
FitDraw_Switch = 0
OriImage_Switch = 1
Data_Total = []
Fit_SaveName = ['NoFit', '']
Data_Cut_Shape_list = ['Rectangle', 'Rectangle']
Data_Cut_Shape = 'RectangleRectangle'
Filter_Size = 2
# %% Modes
File_Mode = 1
Atom_Mode = 0

# %% Initial settings
'''
No need to -1 on row number.
'''
Mat_Data = np.array(range(244, 245 + 1, 1))
# Mat_Data = np.concatenate(([122, 123, 124], np.arange(158, 166 + 1)))
Mat_Data = list(Mat_Data)
# Mat_Data.append(80)
# Mat_Data = np.delete(Mat_Data, np.where(Mat_Data == 5)[0])
# Mat_Data = list(Mat_Data)
Data_Exl = xlrd.open_workbook(Computer_Name + 'Rb87files\\3DSOC_Quench\\Quench.xlsx')
Table = Data_Exl.sheets()[0]
# Matrix_Shift = [0, 0]
clims = [-0.1, 0.5]
FBZ_Len = np.array([104, 52])
# FBZ_Len = 90
Ref_Photo_Num = (20200115, 'line54')
check_photo_number = []
# Wait_Sec = 1
# Photo_Name_Tag='22860test_'
# Photo_Total = ''
# if File_Mode == 0:
#    PhotoName_Pre = 'ascOD'+Photo_Name_Tag
#    PhotoName_End = '.asc'
#    
# elif File_Mode == 1:
#    PhotoName_Pre = Photo_Name_Tag
#    PhotoName_End = 'od.mat'

# %% Calculate all photos

for Row_Num in range(len(Mat_Data)):
    Bad_Pics = []
    Fails = []
    Data_ThisRow = []
    Exl_RowNum = Mat_Data[Row_Num] - 1
    Input_Matrix = Table.row_values(Exl_RowNum)
    Parameters = {}

    # Set parameters
    YearMonthDay = str(int(Input_Matrix[0]))
    Parameters['YearMonthDay'] = YearMonthDay
    Ini_Photo = int(Input_Matrix[1])
    Parameters['Ini_Photo'] = Ini_Photo
    Fin_Photo = int(Input_Matrix[2])
    Parameters['Fin_Photo'] = Fin_Photo
    Repeat = int(Input_Matrix[3])
    Parameters['Repeat'] = Repeat
    V0 = Input_Matrix[4]
    Parameters['V0'] = V0
    Vy = Input_Matrix[5]
    Parameters['Vy'] = Vy
    Omega = Input_Matrix[6]
    Parameters['Omega'] = Omega
    Detuning = Input_Matrix[10]
    Parameters['Detuning'] = round(Detuning, 3)
    DT6_AMP = Input_Matrix[9]
    Parameters['DT6_AMP'] = round(DT6_AMP, 3)
    Cam = Input_Matrix[8]
    Time_Start = Input_Matrix[11]
    Time_Interval = Input_Matrix[12]
    Time_End = Input_Matrix[13]
    if Cam != '22834':
        continue
    Photo_Name_Tag = Cam + 'test_'
    Photo_Total = ''
    if File_Mode == 0:
        PhotoName_Pre = 'ascOD' + Photo_Name_Tag
        PhotoName_End = '.asc'
    elif File_Mode == 1:
        PhotoName_Pre = Photo_Name_Tag
        PhotoName_End = 'od.mat'
    Parameters_str = str(Parameters)
    Parameters_str = Parameters_str.replace('{', '')
    Parameters_str = Parameters_str.replace('}', '')
    Parameters_str = Parameters_str.replace(' ', '')
    Parameters_str = Parameters_str.replace(':', '=')
    Photo = str(Ini_Photo) + 'To' + str(Fin_Photo)
    Photo_Total += Photo + ';'
    # Data matrices
    # if Row_Num == 0:
    #     Time = np.linspace(t_Ini, t_End, (t_End-t_Ini)/Time_Step+1)

    # Addresses
    if File_Mode == 0:
        Folder_Name = r'\\CAMERA-PC\ShareFolder' + '\\' + \
                      YearMonthDay + 'morning\\' + Cam + '\\'
    elif File_Mode == 1:
        Folder_Name = '\\\\calcenter3\\' + 'oriData' + '\\' + YearMonthDay + 'New_Method\\'
    Save_Folder_Name = '\\\\calcenter2\\' + 'oriData' + '\\' + YearMonthDay + 'New_Method\\'
    os.system('mkdir ' + Save_Folder_Name)
    SaveFigure_Address = '\\\\calcenter2\\' + r'Rb87figures\3DSOC_Quench\figures' \
                         + Fit_SaveName[FitBEC_Switch] + '\\' + YearMonthDay + '\\' + Photo + '\\'
    os.system('mkdir ' + SaveFigure_Address)
    SaveData_Address = Computer_Name + 'Rb87Data\\3DSOC_Quench\\' + \
        Fit_SaveName[FitBEC_Switch] + '\\' + YearMonthDay + '\\' + Photo + '\\'
    os.system('mkdir ' + '\"' + SaveData_Address + '\"')
    SaveFigure_Address_All = '\\\\calcenter2\\' + r'Rb87figures\3DSOC_Quench\figures' \
                             + Fit_SaveName[FitBEC_Switch] + '\\' + YearMonthDay + '\\' + 'All_Spin_Polarization' + '\\'
    os.system('mkdir ' + SaveFigure_Address_All)
    SaveData_Address_All = Computer_Name + 'Rb87Data\\3DSOC_Quench\\' + \
        Fit_SaveName[FitBEC_Switch] + '\\' + YearMonthDay + '\\' + 'All_Spin_Polarization' + '\\'
    os.system('mkdir ' + '\"' + SaveData_Address_All + '\"')
    Atom_Center = Quench_Initial_Settings_zdir.AtomCenter_Pos(Date=YearMonthDay)
    Cut_Width = Quench_Initial_Settings_zdir.CutWidths()

    # Calculating
    #    Photo_Ind = Ini_Photo
    SPs_1d = collections.OrderedDict()
    SPs_rened_1d = collections.OrderedDict()
    xs = collections.OrderedDict()
    for Photo_Ind in range(Ini_Photo, Fin_Photo + 1):
        Time = (Photo_Ind - Ini_Photo) * Time_Interval + Time_Start
        plt.close('all')
        Load_Name = Folder_Name + PhotoName_Pre + '%03d' % Photo_Ind + PhotoName_End
        Save_Name = Save_Folder_Name + PhotoName_Pre + '%03d' % Photo_Ind + 'od'
        Fails_All = []
        Data_ThisPic_AtomNum, Data_ThisPic_AtomRegion, ni = \
            bf.Load_Pic_And_Cal(Photo_Ind, Load_Name, File_Mode, Save_Name,
                                SaveFigure_Address + 'All_' + Cam + '_', Atom_Center, Cut_Width, len(Cut_Width) - 1,
                                OriImage_Switch, Rot_Ang=58, Enclose=Data_Cut_Shape_list)
        if Data_ThisPic_AtomNum == 'Error':
            Fails_All.append((Data_ThisPic_AtomRegion, ni))
            Data_ThisPic_AtomNum = 0
            Data_ThisPic_AtomRegion = [np.zeros((Cut_Width[0], Cut_Width[0])), np.zeros((Cut_Width[1], Cut_Width[1]))]

            ni = 0
        Fit_Regions = []
        Fit_AtomNums = []
        TPx = [[], []]
        Sum_0 = np.sum(Data_ThisPic_AtomRegion[0])
        Sum_1 = np.sum(Data_ThisPic_AtomRegion[1])
        Discard_Pic_Threshold = 250
        #        FitBEC_Switch = 0
        #        if Sum_0 < Discard_Pic_Threshold or Sum_1 < Discard_Pic_Threshold:
        #            print('Sum = %.2f and %.2f, no need to fit'%(Sum_0, Sum_1))
        #            FitBEC_Switch = 0
        if FitBEC_Switch == 1:
            for Atom_Ind in range(2):
                AtomRegion = Data_ThisPic_AtomRegion[Atom_Ind].copy()

                plt.imshow(AtomRegion, vmin=-0.1, vmax=0.5, cmap='jet')
                plt.title('%d, %d' % (Photo_Ind, Atom_Ind))
                plt.pause(0.01)
                IsBEC_Threshold = [0.013, 0.020]
                Region_Size = np.shape(AtomRegion)
                AtomCloud_Size = bf.Cal_AtomSize(AtomRegion, [50, 0.20])
                if AtomCloud_Size == 0:
                    AtomCloud_Size = 10000
                Od_Max = np.max(AtomRegion)
                Fi = Od_Max / AtomCloud_Size
                BEC_Judge = bf.Is_BEC(AtomRegion, AtomCloud_Size, IsBEC_Threshold)
                print('BEC_Judge =', BEC_Judge)
                if np.sum(AtomRegion) < Discard_Pic_Threshold and BEC_Judge == 0:
                    print("The atom cloud is too hot. Quit fitting.")
                    Fit_Region = AtomRegion.copy()
                    Fit_Regions.append(Fit_Region)
                    Fit_AtomNums.append(np.sum(Fit_Region) * 112.2)
                    Check_x = np.array(
                        [[1, 0, 0, 100, 100, 0, 3, 0, 30, 30, 0.5], [1, 0, 0, 100, 100, 0, 3, 0, 30, 30, 0.5],
                         [1, 0, 0, 100, 100, 0, 3, 0, 30, 30, 0.5]])
                    Check_x = Check_x.T
                    TPx.append(Check_x)
                    continue
                Fitting_Parameters = np.zeros((11, 3))
                #               0Amp, 1x0, 2y0, 3wx, 4wy, 5rot_ang, 6A2, 7offset, 8wxc, 9wyc, 10Ac
                if BEC_Judge == 1:
                    Od_Lim = 0.3
                    Fitting_Parameters[0, :] = [1, 0, 2]  # Amp
                    Fitting_Parameters[1, :] = [0, -Region_Size[0] / 2, Region_Size[0] / 2]  # x0
                    Fitting_Parameters[2, :] = [0, -Region_Size[1] / 2, Region_Size[1] / 2]  # y0
                    Fitting_Parameters[3, :] = [100, 50, Region_Size[0]]  # wx
                    Fitting_Parameters[4, :] = [100, 50, Region_Size[1]]  # wy
                    Fitting_Parameters[5, :] = [0, -np.pi, np.pi]  # rot_ang
                    Fitting_Parameters[6, :] = [3, 0, 10]  # A2
                    Fitting_Parameters[7, :] = [0, -1, 1]  # offset
                    Fitting_Parameters[8, :] = [10, 0, 20]  # wxc
                    Fitting_Parameters[9, :] = [10, 0, 20]  # wyc
                    Fitting_Parameters[10, :] = [0.5, 0, 1]  # Ac
                    x_Start = Fitting_Parameters[:, 0]
                    x_LowBound = Fitting_Parameters[:, 1]
                    x_UpBound = Fitting_Parameters[:, 2]
                elif BEC_Judge == 0.5:
                    Od_Low = 0.5
                    Od_Lim = np.max((np.max(AtomRegion * 0.5), Od_Low))
                    Size_Start = np.max((AtomCloud_Size, 11))
                    Fitting_Parameters[0, :] = [Od_Lim / 2, 0, Od_Lim]  # Amp
                    Fitting_Parameters[1, :] = [0, -Region_Size[0] / 2, Region_Size[0] / 2]  # x0
                    Fitting_Parameters[2, :] = [0, -Region_Size[1] / 2, Region_Size[1] / 2]  # y0
                    Fitting_Parameters[3, :] = [Size_Start, 10, Region_Size[0]]  # wx
                    Fitting_Parameters[4, :] = [Size_Start, 10, Region_Size[1]]  # wy
                    Fitting_Parameters[5, :] = [0, -np.pi, np.pi]  # rot_ang
                    Fitting_Parameters[6, :] = [Od_Lim / 2, 0, Od_Lim]  # A2
                    Fitting_Parameters[7, :] = [0, -1, 1]  # offset
                    Fitting_Parameters[8, :] = [30, 10, 50]  # wxc
                    Fitting_Parameters[9, :] = [30, 10, 50]  # wyc
                    Fitting_Parameters[10, :] = [Od_Lim / 2, 0, Od_Lim]  # Ac
                    x_Start = Fitting_Parameters[:, 0]
                    x_LowBound = Fitting_Parameters[:, 1]
                    x_UpBound = Fitting_Parameters[:, 2]
                else:
                    Od_Lim = 1
                    Fitting_Parameters[0, :] = [Od_Lim / 2, 0, Od_Lim]  # Amp
                    Fitting_Parameters[1, :] = [0, -Region_Size[0] / 2, Region_Size[0] / 2]  # x0
                    Fitting_Parameters[2, :] = [0, -Region_Size[1] / 2, Region_Size[1] / 2]  # y0
                    Fitting_Parameters[3, :] = [AtomCloud_Size, 0, Region_Size[0]]  # wx
                    Fitting_Parameters[4, :] = [AtomCloud_Size, 0, Region_Size[1]]  # wy
                    Fitting_Parameters[5, :] = [0, -np.pi, np.pi]  # rot_ang
                    Fitting_Parameters[6, :] = [Od_Lim / 2, 0, Od_Lim]  # A2
                    Fitting_Parameters[7, :] = [0, -1, 1]  # offset
                    Fitting_Parameters[8, :] = [7, 0, 40]  # wxc
                    Fitting_Parameters[9, :] = [7, 0, 40]  # wyc
                    Fitting_Parameters[10, :] = [Od_Lim / 2, 0, Od_Lim]  # Ac
                    x_Start = Fitting_Parameters[:, 0]
                    x_LowBound = Fitting_Parameters[:, 1]
                    x_UpBound = Fitting_Parameters[:, 2]
                try:
                    Fit_Region, Check_x = bffit.Main_Gauss(AtomRegion, x_Start, x_LowBound, x_UpBound, Od_Lim,
                                                           'Bimodal')
                except:
                    # Fails.append(('Exl_RowNum and Photo_Ind', Exl_RowNum, Photo_Ind))
                    Fails.append(Photo_Ind)
                    print('Fit failed. Number saved. Check again later.\n')
                    Fit_Region = np.zeros(Region_Size)
                # Fit_Region = bf.Find_Atom_Contour(Fit_Region, 0.05, 1000)
                Fit_Regions.append(Fit_Region)
                Fit_AtomNums.append(np.sum(Fit_Region) * 112.2)
                TPx.append(Check_x)
            Fit_AtomNumsArray = np.array(Fit_AtomNums)
            Data_ThisPic_AtomNum = Fit_AtomNumsArray / np.sum(Fit_AtomNumsArray)
            Data_ThisPic_AtomNum = Data_ThisPic_AtomNum.tolist()
            Data_ThisPic_AtomRegion = Fit_Regions.copy()

        FBZ_Len_1d = FBZ_Len[0]
        Data_Up = Data_ThisPic_AtomRegion[0].copy()
        Data_Dn = Data_ThisPic_AtomRegion[1].copy()
        Data_Up_1d = np.sum(Data_Up, 1)
        Data_Dn_1d = np.sum(Data_Dn, 1)

        Img_Len_Up_1d = Cut_Width[0][1]
        Img_Len_Dn_1d = Cut_Width[1][1]

        DUL_N = 3
        DUL_Range_1d = DUL_N * FBZ_Len_1d
        Data_Up_Large_1d = np.zeros((DUL_Range_1d,))
        Data_Up_Large_1d[(DUL_Range_1d - Img_Len_Up_1d) // 2:(DUL_Range_1d + Img_Len_Up_1d) // 2] = Data_Up_1d
        FBZ_Starts_Up_1d = np.arange(DUL_N) * FBZ_Len_1d

        BZs_Up_1d = []
        # ax.imshow(Data_Up_Large, vmin=-0.1, vmax=0.5, cmap='jet')
        for Starts_Up_1d in FBZ_Starts_Up_1d:
            Ends_Up_1d = Starts_Up_1d + FBZ_Len_1d
            BZs_Up_1d.append(Data_Up_Large_1d[Starts_Up_1d:Ends_Up_1d])
        BZs_Up_1d = np.array(BZs_Up_1d)
        FBZ_Up_1d = np.sum(BZs_Up_1d, 0)
        plt.plot(FBZ_Up_1d)
        plt.title("Time:%dμs Up" % Time)
        plt.show()

        DDL_N = 4
        DDL_Range_1d = DDL_N * FBZ_Len_1d
        Data_Dn_Large_1d = np.zeros((DDL_Range_1d,))
        Data_Dn_Large_1d[(DDL_Range_1d - Img_Len_Dn_1d) // 2:(DDL_Range_1d + Img_Len_Dn_1d) // 2] = Data_Dn_1d
        FBZ_Starts_Dn_1d = np.arange(DDL_N) * FBZ_Len_1d
        BZs_Dn_1d = []

        for Starts_Dn_1d in FBZ_Starts_Dn_1d:
            Ends_Dn_1d = Starts_Dn_1d + FBZ_Len_1d
            BZs_Dn_1d.append(Data_Dn_Large_1d[Starts_Dn_1d:Ends_Dn_1d])
        BZs_Dn_1d = np.array(BZs_Dn_1d)
        FBZ_Dn_1d = np.sum(BZs_Dn_1d, 0)
        plt.plot(FBZ_Dn_1d)
        if np.min(FBZ_Dn_1d) < 0 or np.min(FBZ_Up_1d) < 0:
            Bad_Pics.append(Photo_Ind)
        plt.title("Time:%dμs Down" % Time)
        plt.show()

        if Filter_Size >= 1:
            FBZ_Up_filtered_1d = nd.filters.generic_filter(FBZ_Up_1d, np.nanmean, size=Filter_Size, mode='mirror')
            FBZ_Up_rened_1d = FBZ_Up_filtered_1d[Filter_Size // 2:np.shape(FBZ_Up_filtered_1d)[0]:Filter_Size]
            FBZ_Dn_filtered_1d = nd.filters.generic_filter(FBZ_Dn_1d, np.nanmean, size=Filter_Size, mode='mirror')
            FBZ_Dn_rened_1d = FBZ_Dn_filtered_1d[Filter_Size // 2:np.shape(FBZ_Dn_filtered_1d)[0]:Filter_Size]
            SP_rened_1d = (FBZ_Up_rened_1d - FBZ_Dn_rened_1d) / (FBZ_Up_rened_1d + FBZ_Dn_rened_1d)
            SPs_rened_1d['%03d' % Photo_Ind] = SP_rened_1d
        else:
            SP_rened_1d = 0
            SPs_rened_1d['%03d' % Photo_Ind] = 0
        SP_1d = (FBZ_Up_1d - FBZ_Dn_1d) / (FBZ_Up_1d + FBZ_Dn_1d)
        SPs_1d['%03d' % Photo_Ind] = SP_1d
        xs['%03d' % Photo_Ind] = TPx
        check_photo_number.append('%03d' % Photo_Ind)
        np.savez(SaveData_Address_All + 'Cam_' + Cam + '_photo' + str(
            Photo_Ind) + r'_Spin_Polarization_FBZ' + Data_Cut_Shape + '_1d',
                 SP_1d=SP_1d, SP_rened_1d=SP_rened_1d, Check_xs_up=xs['%03d' % Photo_Ind][0],
                 Check_xs_dn=xs['%03d' % Photo_Ind][1], FBZ_Up_1d=FBZ_Up_1d, FBZ_Dn_1d=FBZ_Dn_1d,
                 Data_Up=Data_Up, Data_Dn=Data_Dn, Data_Up_1d=Data_Up_1d, Data_Dn_1d=Data_Dn_1d,
                 Bad_Pics=Bad_Pics)

        fig, ax = plt.subplots(1)
        plt.plot((FBZ_Up_1d - FBZ_Dn_1d) / (FBZ_Up_1d + FBZ_Dn_1d))
        plt.xlim([0, len((FBZ_Up_1d - FBZ_Dn_1d) / (FBZ_Up_1d + FBZ_Dn_1d))])
        plt.ylim([np.min((FBZ_Up_1d - FBZ_Dn_1d) / (FBZ_Up_1d + FBZ_Dn_1d)),
                  np.max((FBZ_Up_1d - FBZ_Dn_1d) / (FBZ_Up_1d + FBZ_Dn_1d))])

        plt.title('detuning=%.2f, Time:%dμs' % (Parameters['Detuning'], Time))
        plt.show()
        fig.savefig(SaveFigure_Address + 'Cam_' + Cam + '_' + str(
            Photo_Ind) + r'Spin_Polarization_FBZ_color-1to1_det' + Data_Cut_Shape + '_1d' + '.png')
    #        fig.savefig(SaveFigure_Address+'Cam_'+Cam+'_'+str(Photo_Ind)+r'Spin_Polarization_FBZ_color-1to1_det'+Data_Cut_Shape+'.pdf')
    #        Photo_Ind += 1
    plt.close('all')
    #    ii = 1

    #    fig = plt.figure(Row_Num)#, figsize=(20, 1000))
    #    plt.figure(figsize=(20, 10))
    SPs_1d_keys = list(SPs_1d.keys())
    figures = []

    for SPkey_1d_Index in range(len(SPs_1d_keys)):
        if np.mod(SPkey_1d_Index, 20) == 0:
            fig = plt.figure(SPkey_1d_Index // 20, figsize=(16, 13))
            figures.append(fig)
        SPkey_1d = SPs_1d_keys[SPkey_1d_Index]
        Photo_Ind = int(SPkey_1d)
        Time = (Photo_Ind - Ini_Photo) * Time_Interval + Time_Start
        SP_1d = SPs_1d[SPkey_1d]
        #        plt.subplot(11, 15, ii)
        plt.subplot(4, 5, np.mod(SPkey_1d_Index, 20) + 1)
        plt.plot(SP_1d)
        plt.title('%dμs,%03d' % (Time, Photo_Ind), fontsize=16)
        plt.subplots_adjust(top=0.93, bottom=0.019, left=0.019, right=0.973, hspace=0.132, wspace=0.119)
        #        ii += 1
        plt.suptitle('Data_Line=%d,V0=%.3fEr,Vy=%.3fEr,Omega=%.3fEr,Cam=%s,detuning=%.3fEr' % (
                     Mat_Data[Row_Num], V0, Vy, Omega, Cam, Detuning), fontsize=20)

        fig.savefig(SaveFigure_Address_All + 'Data_Line=%d' % Mat_Data[
            Row_Num] + 'Cam_' + Cam + '_' + Photo + 'Spin_Polarization_FBZ_No.%d' % (
                                SPkey_1d_Index // 20) + Data_Cut_Shape + '_1d' + '.png')
    #        fig.savefig(SaveFigure_Address_All+'Cam_'+Cam+'_'+Photo+'Spin_Polarization_FBZ_No.%d'%(SPkey_Index//20)+Data_Cut_Shape+'.pdf')
    plt.show()
    np.savez(SaveData_Address + 'Cam_' + Cam + '_' + Photo + r'Spin_Polarization_FBZ' + Data_Cut_Shape + '_1d',
             SPs_1d=SPs_1d, xs=xs, Fails=Fails, Bad_Pics=Bad_Pics)

    # fig.savefig(SaveFigure_Address_All+'Cam_'+Cam+'_'+Photo+'Row=%d'%Mat_Data[Row_Num]+r'Spin_Polarization_FBZ'+'.png')
    # fig.savefig(SaveFigure_Address_All+'Cam_'+Cam+'_'+Photo+'Row=%d'%Mat_Data[Row_Num]+r'Spin_Polarization_FBZ'+'.pdf')
    np.savez(SaveData_Address_All + 'Cam_' + Cam + '_' + Photo + r'Spin_Polarization_FBZ' + Data_Cut_Shape + '_1d',
             SPs_1d=SPs_1d, xs=xs, Fails=Fails, SPs_rened_1d=SPs_rened_1d, Bad_Pics=Bad_Pics)

# %%
time_end = time.time()
# np.savez('D:\\Rb87figures\\3DSOC_Quench\\figures\\20190430Fails154', Fails=Fails)
print('Time spent: ', time_end - time_start, 's')
time_end_struct = time.localtime(time_end)
print('Start time:', time_start_struct)
print('End time:', time_end_struct)
